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MORBIDITY/MORTALITY  ANALYSIS  OF  COHORT  DATA 


1.  INTRODUCTION 

'5~^,This  report  is  the  culmination  of  investigations  into  modeling  the  natural  history  of  disease  with 
application  to  the  Air  Force  Health  Study_-Tfcree  approaches  to  the  problem  arc  presented  here.  The 
first  approach  extends  the  totally  nonparametric  method  developed  in  Albert,  Gcrtman  and  Louis 
(1978),  Albert,  Gcrtman,  Louis,  and  Liu  (1978)  and  Louis,  Arthur,  and  Hcghinian  (1978).  These  three 
papers  will  be  referred  to  as  Louis  et  al.  The  second  modeling  approach  is  based  on  the  semi- 
parametric  methods  for  analyzing  survival/sacrifice  experiments  presented  in  Dinsc  (1982,  1988)  and 
Portier  and  Dinse  (1987).  This  set  of  papers  will  be  referred  to  as  Dinse  ct  al.  The  third  approach 
combines  the  first  two  modeling  methods  to  obtain  a  semi-parametric  procedure  using  the  disease 
model  of  Louis  et  al.  and  allows  the  inclusion  of  covariate  data. 

In  Section  2  we  present  the  extension  of  Louis  et  al.  to  include  death  as  a  third  time-of- 
occurrcncc.  In  Section  3  we  simplify  the  model  and  obtain  a  full  likelihood  solution.  In  Section  4  we 
show  the  drawbacks  of  attempting  to  extend  the  nonparametric  model  of  Dinse  (1982).  In  Section  5  we 
extend  Dinse  et  al.  to  arrive  at  a  parametric  model.  In  Section  6  we  present  a  partial  likelihood  solution 
to  the  problem  of  determining  the  effect  of  individual  risk  factors  on  death  with  disease  and  death 
without  disease.  In  Section  7  we  present  the  full  likelihood  solution.  In  a  subsection  of  Section  7  we 
also  present  the  formulas  for  a  score  test  to  determine  the  significance  of  risk  factors. 

2.  LOUIS  ET  AL.  EXTENSION  TO  MULTIPLE  EXAMS 

For  multiple  exams,  the  information  available  can  be  written  in  two  parts:  a  basic  Mortal¬ 
ity/Morbidity  vector  denoted  by  MM  and  a  matrix  of  examination  data  denoted  by  E.  Let  T  =  the  age  at 
death  of  a  subject  who  has  died,  X  =  age  at  which  the  preclinical  stage  of  the  disease  begins  and 
Z  =  X  +  Y  be  the  age  at  which  symptoms  appear.  Y  is  the  sojourn  time  in  the  preclinical  stage. 

The  array  MM  is  denoted  by: 

MM  =  (c,DELD,DS,S,C,DR,DC,I), 

where 


e  =  age  of  subject  at  entrance  into  the  study 
DS  =  min(T,  time  of  analysis), 

DELD  =  0  if  DS  =  T, 

=  1  otherwise, 

Z  =  time  of  symptoms  (of  the  disease  of  interest), 

S  =  min(time  to  symptoms,  time  of  analysis,  death  without  symptoms). 
C  =  1  if  subject  is  lost  to  follow-up  at  time  S, 

=  0  if  subject  has  shown  symptoms  al  time  S, 

DR  =  1  if  DS  >  X  and  DELD  =  0, 

=  0  if  DS^  X  and  DELD  =  0, 

DC  =  1  if  DS^X  +  Yand  DELD  =  0, 


=  if  DS  <  X  +  Y  and  DELD  =  0, 


and  I  is  the  number  of  examinations  (or  screens)  the  subject  received.  DR  and  DC  arc  undefined  when 
DELD  =  1.  The  array  E  is  given  by,  E  =  (E(1),...,E(1)),  where  E(i)  =  (U(i),R(i)) ' ,  U(i)  is  the  subject’s 
age  at  the  ith  exam  and  R(i)  =0  if  no  disease  is  present  at  the  ith  exam  and  1  otherwise. 
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The  basic  goal  is  to  estimate  the  joint  probability  density  function  of  (X,Z,T)  in  terms  of  MM  and 
E.  To  illustrate  the  approach,  we  present  here  the  calculation  of  the  likelihood  contribution  of  a  partic¬ 
ular  realization  of  MM  and  E.  Partition  (or  stratify)  (he  positive  real  line  into  intervals  1(1),..., 1(M)  and 
suppose  that  a  subject’s  MM  vector  and  E  matrix  are: 

MM  =  (ee  1(a),  DELD  =  0,  DSe  1(4 ),  Se  I(k),  C  =  1,  DR  =  1,  DC  =  0, 1  =  2) 

and 

UOJelfl,),  R(1)  =  0 
U(2)e(  j2 ),  R(2)  =  1 

MM  =  (a,0,4  ,k,l,  1,0,2) 


Let  us  review  what  information  this  data  represents.  In  Figure  1  we  abuse  notation  and  treat  an  interval 
as  a  point  in  time.  Figure  1  shows  the  disease  and  death  history  of  a  subject  with  this  MM  and  E  data. 


1(a)  I, 

h)  1 

J2) 

4)  I 

:«o 

e  U(l)  U(2)  T  S 


<3 - X - O 


Figure  1. 

Graphical  representation  of  information  available  on  a  subject  with  age  at  entry  (e)  in  the  time 
interval  1(a),  age  at  death  (T)  in  the  interval  1(4 ),  age  at  first  examination  U(l)  in  the  inter  al 
l(jl),  age  at  second  examination  U(2)  in  the  interval  I(j2)>  and  age  at  loss  to  follow-up  (S)  w  the 
interval  I(k). 


Shorthand  totation  would  be 
and 


In  this  special  case,  we  know  Z  > T,  since  C  =  1.  We  also  know  X  <  U(2),  be _ause  R(2)  - 1  (and 
we  assume  there  are  no  false  positives):  however,  if  the  screen  at  U(l)  was  a  false  negative,  X  maybe 
less  that  U(l).  We  also  know  DS  =  T,  the  actual  age  at  death,  since  DELD =0.  Finally,  I  “2  tells  us 
that  this  patient  received  two  exams.  Calculations  for  the  likelihood  contribution  of  this  data  are  (with 
p  =  probability  of  a  false  negative  result): 


h 

P(MM,E)  =  l  P(Xe  I(q),MM,E) 
q  =  l 

M  j2 

=  l  l  P(XeI(q),ZeI(r),MM,E) 

r  =1  +1  q  =  1 

M  j2 

=  £  X  p(PosscreenatIG2)|  «eI(a),X6l(q),ZGl(r),TGl(i),E(l),I  =  2) 

r  =  1  +1  q  =  l 

x  P(e  e  1(a), Xe  I(q),Ze  I(r),Te  1(1  ),E(  1))) 

M  j2 

=  l  l  lx  P(neg.screen at I(jj) |  e 6 1(a), Xe I(q),Ze I(r),Te  1(1  ),I  =  2) 

r=l  +1  q  =  1 

X  P(e  e  1(a), Xe I(q),Ze  I(r),T€  1(1  ),1  =  2) 
x  P(U(l)eI(j1),U(2)eI(j2)  |  «el(a)) 

M  j2 

«  l  l  p w Gl-q)p(e  e  1(a), Xe  I(q),Ze  I(r),Te  1(1  ),I  =  2), 
r  =  l  + 1  q  =  1 

where  u(t)  =  l  t>0,u(t)=0  t^O  and «  indicates  "proportional  to". 
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Note  that  Tel(l )  and  Zel(r),  with  re  {1  4-  is  the  information  given  by 

{  (C  =  l,Se  I(k),Te  1(1  ),DELD  =  0} .  Under  the  assumption  that  examination  times  are  scheduled 
regularly  and  independently  of  (X,Z,T),  as  is  the  case  in  the  Air  Force  Health  Study,  we  have: 

M  J2 

P(MM,E)oc  l  l  pWGl-q)p(a,q,M)P(I  =  2) 

r =1  +1  q  =  l 

M  J2 

«  Z  l  PwGl-q)p(a,q,rye), 
r  =  l  +1  q  =  1 

where  P(a,q,r^ )  =  P(e  e  1(a), Xe  I(q),Ze  I(r),Te  1(1 )).  The  parameters  p(a,q,r,l )  are  the  same  as  those 
of  the  one-examination  model.  Thus,  the  Louis  et  al.  model  extended  to  include  mortality  data  employs 
the  same  parameters  in  the  case  of  two  examinations  as  in  the  case  of  one  examination. 


3.  SIMPLIFIED  NONPARAMETRIC  MODEL 


In  this  simplified  extension  of  Louis  et  al.,  we  reduce  the  model  by  eliminating  X.  We  demon¬ 
strate  the  technique  by  writing  out  the  necessary  equations  in  the  special  case  that  two  exams  have  been 
administered  and  we  use  the  language  of  the  Air  Force  Health  Study,  since  this  form  of  the  model  is 
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specifically  designed  for  ihe  dala  from  (his  study.  Each  subject’s  lime  is  measured  from  his  first  tour  in 
Vietnam.  This  simplified  model  uses  only  the  following: 

T  =  time  of  death, 

Z  =  time  of  symptoms  (of  the  disease  of  interest), 

DS  =  min(T,  time  of  analysis), 

DELD  =  0  if  DS  =  T, 

=  1  otherwise, 

S  =  min(time  to  symptoms,  time  to  analysis  or  death  without  symptoms), 
DELZ  =  1  if  S  =  Z, 

=  0  otherwise, 

I  =  number  of  exams, 

U(j)  =  time  of  the  jth  exam, 

R(j)  =  lifU(j)>S, 

=  0  otherwise. 

The  lime  interval  (0,°°)  is  partitioned  into  intervals  I(1),I(2),...,I(M).  Table  1  lists  the  possible  data  pat¬ 
terns  (indexed  by  rj ). 


TABLE  1.  POSSIBLE  DATA  PATTERNS  IN  THE  REDUCED  LOUIS  ET  AL.  MODEL* 


Pattern  No.  S 

1  k 

2  k 

3  k 

4  k 

5  k 

(. 

7  k 

X  k 

'>  k 

10 

11  k 

12  k 

13  k 

14 

15  k 

10  k 

17  k 

18 


DELZ 

1 

1 

0 

0 

1 

1 

l 

1 

1 

0 

1 

1 

1 

0 

0 

1 

1 


DS 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 

d 


*  Lower  case  letters  (k,  d,  j  j ,  j2  =  1,...,M)  in  the 
in  the  partition  of  (0,  ® ). 


DELD  I 

0  2 

0  2 

0  2 

0  1 

0  1 

0  0 

1  2 

1  2 

1  2 

1  2 

1  2 

1  1 

1  1 

1  1 

1  1 

1  0 

1  0 

1  0 

Z,  DS,  U(l),  U(2) 


U(D  R(D  U(2) 
jl  -  0  \2 

jl  1  J2 

Jl  0  12 

Jl  o 

Jl  1 

Jl  0  J2 

Ji  1  h 

Jl  0  J2 

jl  0  12 

Jl  0  J2 

ii  1 

ji  o 

11  0 

ii  o 


columns  indicate  interval  numbers 


R(2) 

1 

1 

0 


1 

1 

0 

0 

0 


The  goal  of  this  nonparametric  approach  is  to  estimate  the  joint  distribution  of  (Z,  T).  Table  2 
lists  the  likelihood  contribution  of  each  data  pattern  in  terms  of  P(a,b)  =  P(Ze  1(a),  Te  1(b)). 
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TABLE  2.  LIKELIHOOD  CONTRIBUTION  OF  DATA 
PATTERNS  IN  THE  REDUCED  LOUIS  ET  AL.  MODEL 


Data  Pattern  n 
1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

*P(a,b)  =  P(Ze  1(a),  Te  1(b)) 


Likelihood  Contribution 
M 

Z  P(M) 
a  =d+i 


M 

I  P(M) 


l  =d  +  l 

M 

M 

z 

z 

P(q,A) 

q=j2  +  l 

A  =d  + 1 

M 

M 

z 

Z 

P(q>  A. ) 

q=jl  +  l 
M 

A  =d  +  l 

Z  P(M) 

A  =d  + 1 


M  M 

Z  Z  P(q,i) 

q  =  1  A  =d  +  l 
P(k,  d) 

P(k,d) 

P(k,  d) 
d 

Z  P(q,d) 
q=j2+i 

M 

Z  P(q,d) 
q=j2+l 
P(k,  d) 

P(k,  d) 
d 

Z  P(q,d) 

q=jl  +  l 
M 

Z  P(q.  d) 

q  =  d  +  1 
M 

•  Z  P(q.  d) 

q  =  d  +  1 
P(k,  d) 
d  M 

Z  Z  P(q.*) 

q  =  1  A =d  +  1 


We  wish  to  maximize  the  resulting  likelihood  in  terms  of  the  parameters  P(a,  b).  Because  the 
likelihood  function  is  too  complicated  to  maximize  directly,  wc  apply  (as  did  Louis  ct  al.)  the  Estimation 
and  Maximization  (EM)  algorithm  to  obtain  maximum  likelihood  estimates  of  the  P(a,  b).  To  invoke 
the  EM  algorithm  we  need  to  compute  the  conditional  probabilities  P(Ze  1(c),  Te  1(0  |  V ),  the  proba- 
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bility  that  Ze  1(c)  and  Te  1(0  given  that  the  subject  has  displayed  data-pattern  tj  .  Using  the  definition 
of  conditional  probability  we  see  that 

P(Ze  1(c),  Te  1(f)  |»?)  =  P(Ze  1(c),  Te  l(0,t?  )/P(ij ),»?  =1,2,...,18, 

where  P(»7 )  is  the  likelihood  contribution  of  the  data  pattern  T) .  To  compute  these  conditional  proba¬ 
bilities  we  first  compute  the  probabilities  of  a  subject  having  Ze  1(c)  (i^symptoms  occur  in  the  interval 
1(c))  and  Te  1(0  (death  occurs  in  interval  1(0)  and  the  subject  displays  data  pattern  tj  .  Let 

5(t)  =  1  if  t  =  0 
=  0  otherwise, 

and  u  (t)  =  1  if  t  >  0 

=  0  otherwise 

Table  3  gives  these  probabilities  for  each  t?  ,t?  =  1,2,...,  18. 


TABLE  3.  THE  PROBABILITIES  THAT  SUBJECTS  SYMPTOMS  APPEAR  AT  TIME  S  IN 
THE  INTERVAL  1(C)  AND  THAT  THE  SUBJECT’S  TIME  OF  DEATH  (T)  IS  IN  THE 
INTERVAL  1(F),  GIVEN  THAT  THE  SUBJECT  HAS  DATA  PATTERN  rj 


tL  (S.  DELZ,  PS.  DELD.  I,  U(l).  Rfl).  U(2),  R(2)) 


1  (k,  1,  d,  0, 2,  jj,  1,  j2, 0) 

2  (k,l,d,0,2,j1,l,j2,l) 

3  (k,0,d,0,2,j1,0,j2,0) 

4  (k,  0,  d,  0, 1,  j  j,  0,  -,  -) 

5  (k,  1,  d,  0, 1,  jx,  1,  -,  -) 

6  (-,  -,  d,  0, 0,  -,  -,  -,  -) 

7  (k,  1,  d,  1, 2,  jj,  0,  j2, 1) 

8  (k,  1,  d,  1, 2,  jj,  1,  j2, 1) 

9  (k,  1,  d,  1, 2,  jj,  0,  j2, 0) 

10  (-,  1,  d,  1,  2,  jj,  0,  j2, 0) 

11  (j2,0,d,  l»2,ji,0,j2,0) 

12  (k,  1,  d,  1, 1,  jj,  1,  -,  -) 

13  (k,  1,  d,  1, 1,  jj,  0,  -,  -) 

14  (-,  1,  d.  1, 1,  jj_,  0,  -,  -) 

15  (d,  0,  d,  1, 1,  jj,  0,  -,  -) 

16  (k,0,d,  1,0, 

17  (k,  1,  d,  1,0, 

18  (k,l,d,  1,0, 


P(Ze  1(c),  Te  KO.  n) 
P(c,0w(f-d)5(k-c) 
P(c,f)w(f-d)<5  (k-c) 

P(c,f)w  (c-j2)w  (f-dl 
P(c,0w(c-j1)w(f-d) 
P(c,f)w(f-d)5  (k-c) 

P(c,f)<J  (f-d) 

P(c,05  (f-d)5  (k-c) 

P(c,f)5  (f-d)5  (k-c) 

P(c,05  (f-d)5  (k-c) 
P(c,0w(c-j2)5  (f-d)(l-«(d-c)) 
P(c,0«(c-d)5  (f-d) 

P(c,05  (k-c)5  (f-d) 

P(c,0<5  (k-c)5  (f-d) 
P(c,f)w(c-j1)5  (f-d)(l-w(d-c)) 
P(c,0«(c-d)5  (f-d) 
P(c,f)w(c-k)5  (f-d) 

P(c,05  (k-c)6  (f-d) 
P(c,0(l-^(d-c))5(f-d) 


The  maximum  likelihood  estimates  of  the  P(a,  b)  are  obtained  by  using  the  following  EM  algo¬ 
rithm.  Let  n(f? )  =  the  number  of  subjects  with  data  pattern  rj ,  T]  =  1,2..., 18. 

Step  1:  Choose  a  set  of  initial  probabilities: 

P(°)  (a,b),  a,b  =  1,...,M. 

Step  2:  Given  the  vth  estimates  p(v)(a,b),  v  =0,1,2,...,  compute 
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nC1'  +  1)(a,b)  —  £  n(t7)p(y)(a(b|  7),  a,b  =  1,...,M, 

V  =1 

where  p(v)(a,b|  7)  is  obtained  by  using  p(v')(a,b)  in  the  formulas  of  Tables  1  and  3  and 
the  formula  for  P(Z f.  1(a), Te  1(b)  |  7 ). 

Step  3:  Compute  p(v  +  *)(a,b)  =  n(v  +  *)/N,  a,b  =  l...,M,  where  N = £  n(7  ). 

V 

Step  4:  Repeat  Steps  2  and  3  until  the  p(v  )(a,b)  converge. 


The  resulting  limits  are  the  maximum  likelihood  estimates  of  the  P(a,b). 


4.  TWO-EXAMINATION  EXTENSION  OF  DINSE  (1982)  AND  ITS  DRAWBACKS 


Dinse  (1982)  provides  a  method  for  analyzing  data  composed  of  examination  and  death  time 
information  when  subjects  are  examined  at  most  one  time.  Below  we  provide  an  example  of  an  exten¬ 
sion  of  this  work  that  allows  more  than  one  examination.  This  development  is  exactly  that  of  Dinse 
(1982)  except  that  the  number  of  examinations  is  2  rather  than  1.  This  approach  illustrates  the  fact  that, 
unlike  the  extension  of  Louis  et  al.  in  Section  2,  increasing  the  number  of  examinations  by  one  increases 
the  number  of  parameters  in  the  model  by  a  number  equal  to  the  product  of  the  number  of  death  times 
and  the  number  of  disease  states.  This  increase  in  parameters  makes  the  method  too  cumbersome  for  a 
study  like  the  Air  Force  Health  Study  in  which  there  are  to  be  up  to  six  examinations.  The  full  data  case 
is  sufficient  to  demonstrate  the  drawback  of  attempting  to  use  a  two-examination  extension  of  Dinse. 

Define  the  following  notation: 


T  =  age  at  death, 

B(T)  =  disease  state  at  death, 

E(i)  =  age  at  the  ith  exam,  i  =  1,2, 

B(i)  =  disease  state  at  the  ith  exam. 

The  full  data  case  involves  no  incomplete  pairs  among  (T,B(T)),  E(i),B(i)),  although  one  or 
more  of  (E(i),B(i))  may  be  missing.  By  design  E(l)  <  E(2),  so  that  there  are  three  possibilities  for  an 
individual:  T  <  E(l)  <  E(2),  E(l)  <  T  <  E(2),  or  E(l)  <  E(2)  <  T.  If  T  <  E(i),  then  the  ith  examination 
was  not  performed  (of  course). 


If  we  follow  the  nonparametric  approach  of  Dinse,  we  would  take  each  possible  event  and  write 
out  its  likelihood  contribution  as  a  product  of  conditional  and  unconditional  probabilities.  These  condi¬ 
tional  and  unconditional  probabilities  would  then  be  considered  parameters  to  be  estimated.  The 
infeasibility  of  this  approach  for  the  Air  Force  Health  Study  can  be  best  illustrated  by  writing  out  the 
likelihood  contribution  for  an  individual  having  had  two  exams  (ie,E(l)  <  E(2)  <  T).  Using  the  time 
intervals  1(1)  <  1(2)  <  •  •  •  <  I(M)  of  the  Louis  et  al.  extension,  the  likelihood  contribution  would  be: 

P(Tg  1(a), B(T)  =  b(T),E(l)e  1(c), B(l)  =  b(l),E(2)e  1(d), B(2)  =  b(2)) 

=  P(B(T)  =  b(T)  j  Te  1(a), E(l)e  1(c), B(l)  =  b(l),E(2)  =  I(d),B(2)  =  b(2)) 
x  P(Te  1(a)  |  E(l)e  1(c), B(l)  =  b(l),E(2)  =  I(d),B(2)  =  b(2))  x  P(B(2)  =  b(2)  |  E(l)e  1(c), 
B(l)  =  b(l),E(2)  =  1(d))  xP(B(l)  =  b(l)  |  E(l)€  I(c))P(E(l)e  I(c),E(2)e  1(d)). 
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The  first  factor  requires  a  parameter  for  every  combination  of  death  time,  exam  times  and  disease 
states.  Essentially,  this  parameter  forces  the  data  to  be  stratified  according  to  each  such  combination. 

This  problem  illuminates  the  difference  between  survival/sacrifice  experiments  and  epidemiolog¬ 
ical  studies.  In  a  survival/sacrifice  experiment  the  sacrifice  is  the  first  and  only  examination  and  coin¬ 
cides  with  the  ending  of  information  on  that  animal.  However,  in  an  epidemiologic  study  with  repeated 
exams,  the  subjects  may  continue  to  live  after  the  first  and  second  exams  and  their  disease  history  will 
take  one  of  many  paths.  To  make  inferences  about  these  paths,  we  must  have  a  sufficient  number  of 
both  cases  and  controls  taking  each  path.  A  large  number  of  paths,  such  as  would  be  the  situation  for 
two  or  more  exams,  would  require  a  very  large  and  therefore  infeasible  number  of  cases  and  controls. 
Clearly,  this  is  an  unrealistic  approach  for  the  Air  Force  Health  Study.  An  alternative  and  feasible 
approach  is  to  invoke  a  parametric  model,  described  in  the  next  section. 

5.  A  PARAMETRIC  MODEL 


Now  let  X  denote  the  time  to  the  first  event,  either  onset  of  the  disease  of  interest  or  death 
without  the  disease  of  interest.  Let  T  denote  the  time  to  natural  death.  Define 

Y(t)  =  1  if  disease  is  present  at  time  t, 

=  0  otherwise. 

Suppose  there  are  J  distinct  death  and  examination  times.  Let  V  be  a  vector  of  covariates 
measured  on  each  subject.  Define 

aj(  V)  =  number  of  subjects  who  died  with  the  disease  at  tj  and  had  covariate  vector  V. 

bj(V)  =  number  of  subjects  who  died  without  the  disease  at  tj  and  had  covariate  vector  V. 

mj(V)  =  number  of  subjects  examined  at  tj  who  were  disease  free  and  had  covariate  vector  V. 

nj(V)  =  number  of  subjects  examined  at  tj  who  had  the  disease  and  had  covariate  vector  V. 

Define  the  following  hazard  functions: 

Xv(t)  =  lim  P(t£X<t  +  e,Y(X)  =  l|  X£t,V  =  v)/e, 
e-0 

/?v(t)  =  lim  P(t2X<t  +  e,Y(X)  =  0|  Xi>t,V  =  v)/e, 
e-0 

=  lim  P(t^T <t  +  e  j  T£ t, Y(t)  =  0, V  =  v)/e , 
e-0 

av(t)  =  lim  P(t£T<t  +  e  |  T£t,Y(t)  =  l,V  =  v)/e. 
e-0 

These  functions  extend  Dinse’s  (1988)  expressions  to  include  covariates.  It  follows  directly  from 
Dinsc  (1988)  that  X  v(t),  £v(t)  and  av(t)  can  all  be  written  in  terms  of  the  following  three  functions: 


8 


and 


h^t)  =  lim  P(t^T<t  +e  ]  T£t,V  =  v)/e, 
e-0 

71  v(0  =  P( Y(t)  =  1 1  t, V  =  v), 


Pv(0  =  P(Y(t)  =  1 1  T  =  t,V  =  v). 

Following  Dinse  (1988),  we  define  the  lethality  function  as  ry(t)  =  ay(t  )/£y(t).  It  follows  from  Dinse’s 
(1988)  expression  (4)  that  rv(t)  =  (pv(t)/(l-pv(t)))/('7tv(t)/l-itv(t))). 


We  now  introduce  the  following  model:  Let  £ ,  v  and  7  be  vectors  of  coefficients  and  £  g(t), 
vg(t)  and  7  g(t)  be  functions  of  time.  We  model  hj(t)  =  hg(t)e£  v,  a  proportional  hazards  model,  and 
7tv(t)  and  pv(t)  with  logistic  models, 


and 


ev0(0  +  v  'v 
l  +  e^flW+v  'v 


Pv(0  = 


670(0  + 7  'v 
l  +  e7Q(0  +  7  'v 


It  follows  that  the  survival  function  Sy(t)  =  Py(T  >  t)  can  be  written  as  a  power  of  a  baseline  survival  function 
determined  by  hg(t), 


We  also  see  that 


Sv(t)  =  exp{-/o  h^ujdu} 

=  (exp{-/  hg(u)du}) 
0 


e£'v 


7  n(0*vft(0  +  (7  '-v  ')v 


r  (t)  =  e 1  0W  0 


so  that  the  ratio  of  lethality  for  covariate  vectors  V  =  vj  to  V = V2  is 


rvl(0/rv2(0  =  e 


(7  -V  ')(vrv 2) 


The  theory  in  section  3  of  Dinse  (1988)  together  with  the  above  models  imply  that  the  logiikelihood  is 
the  sum  over  all  vectors  z  of  the  following  expressions. 


9 


J  tj 

Ll(hv)  =  I  {(aj(v)  +  b:(v))(logho(tj)  +  £ 'v)-(aj(v)+b:(v)  +  mj(v)+iij(v))e«  VJ  hQ(u)du} 
j=l  1  1  0 

l^Oty)  =  I  {nj(v)(i/0(tj)+i'  'v)-(nj(v)  +  mj(v))log(l  +  ev0^tj^+l'  Z)} 

j  =  l 

L>3(Pv)  =  I  {bj(v)(70(tj)  +  7  'v)-(aj(v)+bj(v))log(l  +  e'r0^ti^+7  Z)}. 

j  =  l 

The  resulting  loglikelihood  is  L  =  £  [Lj(hv)  +  L2(iHv)  +  L3(pv)]. 


The  baseline  information  hg(t),  vq(0  and  7g(t)  will  be  known  if  the  exposed  group  is  being 
compared  with  a  population  (Breslow,  Lubin,  Marek,  and  Langholz  (1983)).  In  a  two-sample  study, 
however,  this  baseline  information  will  not  be  known  and  thus  must  be  estimated.  Such  estimates  are 
derived  via  modeling.  Survival  times  are  frequently  assumed  to  follow  a  Weibull  model  (Dinse  1988). 
Dinse  and  Lagakos  (1982)  suggest  modeling  the  logit  of  the  prevalence  rate  by  a  low-order  polynomial 
in  time.  Taking  these  suggestions,  we  model 

(A  i-l 

hg(t)=^oMlf  . 


v0(t)=v0  +  vxt+v2t2 

and 

Tf0(t)  =  80  +  8lt  +  82t2- 

The  resulting  parameters  are  estimated  by  maximizing  the  full  likelihood. 


5.1  Some  Special  Cases 


In  the  following  special  cases  we  assume  v  is  a  group  indicator,  v  =  1  for  exposed  and  0  for 
controls. 

Case  1:  pv(t)  =  p0(t),  7t  v(t)  =  7tg(t)  and  hg(t)  is  known. 


This  is  the  case  assumed  in  Breslow,  Lubin,  Marek,  and  Langholz  (1983).  In  this  situation  the 
only  difference  between  exposed  and  controls  is  in  survival.  The  loglikelihood  in  this  situation  is 
L  =  £vLj(hv).  The  hypothesis  of  interest  is  Hq:  £  =0.  The  corresponding  score  test  (Rao,  1973) 
which  is  based  on  the  statistic  S  =  (d  Ud £  |  £  =g)^/(-3^U3^^|  £  _g).  In  this  case, 


>L/3£  |  £  =0  = . I  {(a;(l)  +  b:(l))-(a:(l)  +  b:(l)  +  n:(l)) /\(u)du}, 
j  =  l  0 


=  O-E 
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where  0=  £  (a:(l) +  b:(l))  =  number  of  deaths  among  the  exposed  group  and 

j  =  l 

J  P 

E  = .  I .{(aj(l)  +  b;(l)  +  m;(l)  +  n:(l))  f  Jhg(u)du} 

j  =  1  J  1  1  J  0 

=  expected  number  of  deaths  among  the  exposed  group 
under  the  control  hazard, 

and  (-3  2 Ud  £2)  |  c  =y  =  E.  Rao’s  score  statistic  for  testing  Hg:  £  =  0  is  thus  (0-E)2/E.  This  is  the 
same  test  statistic  obtained  by  Breslow,  Lubin,  Marek,  and  Langholz  (1983)  in  a  similar  situation. 
Case  2:  £  =  0, 7  =  0  with  v  g(t)  and  7  g(t)  known. 


This  is  the  case  in  which  exposed  and  controls  differ  only  in  the  prevalence  of  disease  among  the 
living,  with  the  prevalence  among  the  living  and  prevalence  at  death  assumed  known  for  the  controls.  In 
this  case  the  observed  number  of  diseased  subjects  among  the  living  exposed  is: 

J 

0  =  1  n:(l), 

j  =  l 

while  the  expected  number  of  disease  subjects  among  the  exposed  under  the  control  prevalence  rate  is: 

E-  1  {mj(  1)  +  nj(l))(e*-  >  +  "/(!  +  ^0  C  j )  +  ^ )). 

i-1 

Hence,  (3  Udv  \  v  =0)  =  OE  and  (32L/3i/^|  ^  =q)  =  -E.  Thus  the  score  statistic  for  testing  Hq: 
v  =0  is  again  of  the  form  (0-E)2/E. 

Case  3:  £  =  0,  v  =  0  with  v  g(t)  and  7  g(t)  known. 


In  this  case  the  only  difference  between  exposed  and  controls  is  the  prevalence  of  disease  at 
death,  with  the  control  prevalence  known.  It  again  follows  that  the  score  statistic  for  testing  Hq:  7  =  0  is 
(0-E)2/E,  where 


and 


J 

O  =  £  b:(l)  =  observed  number  of  diseased  cases  among  those  dying  of  natural  causes, 

j  =  l 


E  =  l  (ajW  +  bjOJXe'1 
j  =  l 


)). 


5.2  The  General  Case 


In  the  general  case  that  the  baseline  functions  cannot  be  assumed  known  and  there  may  be 
covariate  effects  for  survival  and  both  prevalences,  we  assume  the  baseline  models  in  (1).  The  loglikeli- 
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hood  is  the  sum  over  z  of  Lj(hz),  L^Oty)  and  (l-3(pv)  with  (1)  substituted  for  h^t)  and  i/q(0  and  7q(0- 
The  resulting  likelihood  is  then  maximized  to  obtain  the  parameter  estimates. 

6.  A  PARTIAL  LIKELIHOOD  MODEL 


Again,  let  Z  indicate  the  time  symptoms  occurred  and  let  T  indicate  the  time  of  death.  In  this 
section,  we  extend  the  Dinse  (1982)  model  in  terms  of  Z  and  T  and  then  use  this  model  to  build  a  likeli¬ 
hood.  The  events  for  which  we  need  likelihood  contributions,  numbered  I,  II,  III,  IV,  and  V  as  in  Dinse 
(1982),  are  all  the  combinations  of  censoring  situations  on  Z  and  T.  Table  4  lists  these  events  and  their 
likelihood  contributions. 


TABLE  4.  OBSERVABLE  EVENTS  AND  THEIR  LIKELIHOOD  CONTRIBUTIONS  FOR 

THE  DINSE  (1982)  EXTENSION 


Event  Tvoe 

Event 

Likelihood  Contribi 

I 

Z  >  t,T  >  t 

P(T  >  1 1  Z  >  t)P(Z  >  t) 

II 

Z<t,T>t 

P(T  >  1 1  Z  <  t)P(Z  <  t) 

'd 

III 

Z  >  t,T  =  t 

—  P(T>t)  Z>t)P(Z>t) 

-d 

IV 

Z  <t,T  =  t 

—  P(T>t|Z<t)P(Z<t) 

V 

T>t 

Sum  of  likelihood  contributions  I,  II 

Notice  that  we  are  assuming  that  the  actual  time  of  occurrence  of  symptoms  is  never  known.  The 
model  can  be  easily  extended  to  include  such  cases  if  necessary.  Notice  also  that  the  events  depict 
knowledge  at  a  point  in  time  T  =  t.  This  time  T  would  be  the  time  of  the  last  available  knowledge  about 
the  subject.  Thus  T  can  be  the  time  of  death,  the  time  of  loss  to  follow-up,  or  time  of  analysis. 


Let  V  indicate  a  vector  of  covariates.  We  define  the  following  conditional  hazards. 

A  (t  |  Z< t,V  =  v)  =  lim  P(t£T<t+e  j  Z<t,V  =  v,T£t)/e 
e-»0 


A  (t  |  Z>t,V  =  v)  =  lim  P(t£T<t+e  [  Z>t,V  =  v,T£t)/e. 
e-0 


We  also  define  the  conditional  baseline  hazards  as: 
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X1(t)  =  X(t|Z<t,V=0) 

and 

X2(t)  =  X(t|Z>t,V  =  0). 
Finally,  we  model  the  conditional  hazards  of  T  as: 

X(t|  Z<t,V  =  v)  =  X1(t)e/9  'v 


and 

Mt  |  Z  >  t,V  =  v)  =  X  2(t)ea  v, 

where  0  and  a  are  vectors  of  parameters  which  reflect  the  effect  of  the  covariates  in  V  on  the  hazard  of 
death  for  a  person  with  and  without  the  disease,  respectively. 

Let  t(]j  <  t(2)  <  •  •  •  <  t^  be  the  ordered  death  times  of  those  dying  with  the  disease.  Denote  by 

Vq  the  vector  of  covariates  for  the  subject  dying  at  tQ.  Similarly,  let  t^  <  t(2)  <  •  •  •  <  t^  be  ordered 

death  times  of  those  dying  without  the  disease.  Let  Vq  be  the  vector  of  covariates  for  the  subject  dying 

at  tQ.  Define  Rq  ^  and  R^;  ^  to  be  the  set  of  indices  for  subjects  at  risk  at  Iq  j  and  t^j  respectively, 

jj  =  l,2,...,k,  j2=  l,2,...,k.  Standard  partial  likelihood  arguments  (Lawless  (1982),  pp356-357)  lead  to  the 


partial  likelihood 


Uptj)  =  ( 


k  e  ^  v0)  ^  ea '  v(z) 

il“l  l  eW  (  fell  l  e“V' 


XeR(Jl)  iG^(j2> 


(2) 


This  partial  likelihood  is  the  product  of  two  standard  partial  likelihood  functions.  Thus,  it  can  be 
maximized  using  standard  methods  on  each  factor  to  derive  estimates  of  a  and  0. 


7.  A  FULL  LIKELIHOOD  MODEL 


The  partial  likelihood  solution  to  (2)  allows  us  to  compare  the  effect  on  death  rates  of  a  covariate 
for  subjects  with  and  without  the  disease.  However,  there  are  other  descriptors  for  the  disease  which 
we  cannot  estimate  without  further  modeling.  For  example,  the  lethality  of  the  disease  at  age  t  for 
subjects  with  covariates  V  =  v  can  be  defined  as: 

X  (t  |  S<t,  V =v) 

"  X  (t  |  S  >  t,  V  =  v)  ■ 

In  terms  of  the  model  in  Section  2  we  have: 


>2(0 

so  that,  without  a  trivializing  assumption  on  the  baseline  hazards,  we  can  make  no  inference  about  rv(t). 
The  solution  to  this  problem  is  to  model  the  baseline  hazards.  In  this  section,  we  give  a  full  likelihood 
solution  to  the  discretized  form  of  this  problem.  This  solution  will  be  useful  in  analyzing  grouped  data. 
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Let  1 2,t2,-..,t|c  be  the  set  of  possible  death  times.  Define: 

Xlj  =  P(T  =  tj|'n>tj,S£tj),j  =  l,...,k 

and 

X  2j  =  P(T  =  tj|T^  tj,S  >  tj),  j  = 

Let 

Dj(d)  =  The  set  of  labels  of  individuals  dying  at  time  tj  and  with  the  disease  (S£  tj) 

Cj(d)  =  The  set  of  labels  of  those  with  the  disease  censored  at  tj. 

Dj(3)  =  The  set  of  labels  of  subjects  without  the  disease  who  died  at  tj. 

Cj(d)  =  The  set  of  labels  of  subjects  without  the  disease  who  were  censored  at  tj. 

Cj  =  The  set  of  labels  of  persons  whose  death  times  were  censored  at  tj  and  for  whom  we 
have  no  disease  information. 


(3) 


(4) 


If  there  are  no  censored  death  times  without  disease  information,  then  L(X  ,ajt)  is  the  product 
of  two  standard  likelihoods  for  the  discrete  proportional  hazards  model.  Substituting  (3)  into  (4)  and 
taking  the  logarithm,  allows  us  to  write  the  log  likelihood  as: 

k 

logL(7,a^)  =.1^1  log(l-exp(-exp(7 Xj  +  0  'vj)))-  £  cxp(7lj  +  yS  'v^) 

J  ieDj(d)  ieCj(d) 


Following  the  technique  of  Kalbfleisch  and  Prentice  (1980),  pp  98-102,  we  model: 

7  ij  ~  h>g((-log(l-X  jj)) 

and 

72j  =  log((-log(l-X2j)). 

Using  the  discrete  form  of  the  proportional  hazards  model,  we  have: 


The  resulting  likelihood  is: 


P(T  =  tj,S<:  tj  1  T*  tj,V  =  v)  =  1-(1-X  jj)e 
P(T  =  tj,S  >  tj  |  T£  tj, V  =  v)  =  1-(1-(1-X  2j))e° 


k 

L(x,a^)=n  {  n  [i-(i-Xjj)c^  'vii  n  (i-x^  'v*  n  [i-(i-x2j)ea'vi! 

J_1  ieD;(d)  XeC;(d)  ieD:® 


xn  (i:x2j)eQt  n  [i-(i-x  !j)e  ^  v-4  +(i-x2j)**'v*]}, 


^eCj(d) 


where  X  (X  ii>***>^  ik»^  21’****^  2k^ 
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+  X  log(l-exp(-exp(72j  +  a 'Vj)))-  X  exp(72j +a'vje) 

ieDj(d)  XeCj(d) 


+  X  log(exp(-exp(7  ij  +  0  '  v* ))  +  exp(-exp(7  2j  +  a  '  ))) } . 

AeCj 

To  form  the  score  test  for  H q.0  =0,a  =0,  we  compute  for  j  =  l,...,k,  jj  =  l,2,...,k,  j2  =  1,2,..., k,  v  =  1,2  and 
r,q  =  l,...,p, 

[5  r  j = ^ — iogL}, 

*7ij 

[a^]  =  {— |-iogL}, 

°P  r 


I  <9  <*]  =  {- - IogL}, 

daT 


~  -3  2logL 

[-3  2/9]  =  { - — 

dfirWq 


-32|ogL 

daTdaq 


^  -o 

[-8  2a]  =  {-  - - }, 


-  -3  2logL 

[-3  2/Sa]  =  {———}, 

3ar3^q 

~  -3  2logL 

°7ij3pr 

-  -3  2logL 

[-a2ria]-{- — — — }, 

0  7  ij3  ar 

and 

-  -3  2logL 

[-a  2r  ,r  2|  -  f- — }. 

aTV,% 

[3  T ;]  is  a  k  x  1  vector  with  jth  component - IogL,  [-3  ^0ot]  is  a  p  x  p  matrix  and  [-3  2r  is  a  k  x  p 

5  7 :: 

~  ~  *  3  IogL 

matrix,  i  =  1,2.  Define  T  to  be  the  vector  of  7 the  solutions  to - =  0  where  the  derivative  is 


evaluated  at  a  =0  and  0=0.  The  components  of  T  can  be  written  simply  with  more  notation. 
Define: 
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and 


ajj  =  The  number  of  indices  in  Dj(d), 
a2j  =  The  number  of  indices  in  Dj(d), 
bjj  =  The  number  of  indices  in  Cj(d), 
t>2j  =  The  number  of  indices  in  Cj(d), 

cj  =  The  number  of  indices  in  Cj. 


Finally  define: 

pvj  =  exp[-exp(71,j)J,  j  =  v  =  1,2. 


With  this  notation  the  components  of  T  are  the  solutions  to  the  following  k  pairs  of  quadratic  equations. 


and 


j  =  l,...,k.  Define: 


Plj(alj  +  blj +  cj)  +  Plj(-blj'cj)  +  PljP2j(alj +  blj)P2jblj  =  0 
P2j(a2j  +b2j  +  cj)  +  P2j(-b2j-cj)  +  PljP2j(a2j +  ^p-Plj^j  =  °> 


and  Y.  as  the  2(k  +  p)x2(k  +  p)  symmetric  matrix 


[-a2rij 

l-32r1r2j 

H2iVl 

[-3  2 I>] 

( -^  2r  2J 

l-a2r^i 

[-a2??*] 

[-a  2/3] 

[-a2P  a] 

[-32a] 

2 

Then,  the  score  test  for  Hq:  fi  =  0,  a  =  0  is  the  X2p  statistic 

X2  =  (d  L)o£o(d  L)o, 

A 

where  the  subscript  "0"  indicates  evaluation  at  (r,ce^)  =(r,0,0). 
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7.1  Computational  Formulas  For  The  Score  Test 

The  following  formulas  are  all  evaluated  at  o  =0, 0  =0.  The  rth  component  of  the  vector  of 
covariates  for  the  ith  subject  is  denoted  by  Z^. 

j-1 . K,  v  =  1,2. 

^  7  v  i  1‘Py  j  Pij +  P2j 

~~~  =  I  {  Pl^°gPlj-  l  Zir  +  logpij  l  ziT+  PlL(^gPll  £  zlr},  r  =  l . p. 

3/9r  j-1  1-Pij  ig  Cj(d)  1 6  Cj(d)  plj  +  p2jleCj 

aiogL^  ^log^  z  ^  £  Z,r+^E2L  X  Zir},  r  =  1 . p. 

aar  J  =  1  1-P2J  feD.^  ieCj(a)  Plj  +  p2j  *eCj 


-a2logL  Plj5ogp1j  +  (logp1j)p1j-p:ijlogp1; 


2  ~  -°lj(' 

^lj 


(i-Pij) 


bijlogpij 


.  _  ,  Pl^ogPij  +  Pi^ogPli  +  PliPZiOoKPli)  x 

+  Cj(  “  )> 

(Plj  +  p2j> 

-32logL 

and  for  j  =  - -  0,  ')*  s. 

a7ljafls 

[-3  2r  j]  is  a  diagonal  matrix.  Similarly,  [-3  2T 2]  is  a  diagonal  matrix  whose  jth  diagonal  element  is: 


<' W  P2j-i4i'°gP2i ,  ,  , 

2.  a2j*  n.^2  )  +  P2j,0gP2j 

7  2] 


(1-P2j)^ 


P2ilogP2j  +  PliP2ilogP2j  +  PliP2i(logP2j) 
(plj  +  p2j)2 


for  j  =  l,.,.,k.  Additionally, 


Z  zir 

3^r37ij  (1-Plj)  i£ D:(d)  J  leC:(d) 


t  PlilogPli(Pli  +  P2j-P2ilogP2i)  y 

<Plj  +  P2j)2  4eC:  lr 
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-3  ^logL 
d  &  rd  y  2j 


-3  2logL 

aara71j 


PliP2iloSP2ilogPlj 
(Plj  +  P2j)2 


■  l  Hr 

AeCj(d) 


-3  2logL 

d0rd7  2j 


j~rL-=  Pl^ 7^1  *****  l  Hr 

3ard7ij  (1-P2j)  ieDjO)  ieCj(a) 


.  P2ilo8P2i(Pli  +  P2i-Plilo«Pli)  r 
+  /  ,  s2  L  Hr 

(Plj  +  P2j>Z  ieQ 


J 


-3  2IogL  *  ,  PijlogPij(l  +  IogPii-Pii) 


=  y  { 

3/Wq  jSi 


tt-pljr 


I  Zir^iq 


feDj(d) 


Plil°gPli(Pli  +  P2i  +  P2i1°gPli)  r  _ 

-‘ogPij  l  *lqHr - 1 - ‘  \2  -  E 

-  -  -  -  •'  Cpij +  P2j)  <•  -  ~ 


jeeCj(d) 


ieC: 


-3  2logL 
3  cer3  ctq 


V  ,  P2il°gP2i(1-P2i  +  1°gP2i)  v  7. 

"j«  <H& 


~logP2j  I  ZJ?  r  zi  q  + 

XeCj(d) 


P2jlogP2j(Pli +  P2j  +  PljlogP2j) 
(Plj  +  P2j)2 


I 

£eCj 


H  rzlq}- 


-3  2logL 
3ar3^q 


I’1  (Plj  +  P2j>  ieCj 


8.  DISCUSSION 


There  remains  much  work  to  be  done.  These  three  models  must  be  compared.  The  most  useful 
model  seems  to  be  the  semiparametric  partial  likelihood  model,  because  software  already  exists  for  the 
estimation  of  its  parameters.  However,  the  simplified  extension  of  Louis  et  al.  settles  down  somewhat  as 
the  number  of  screenings  increases.  This  suggests  that,  if  sample  sizes  are  large  enough  for  one  screen, 
they  are  also  large  enough  for  multiple  screens.  Of  course,  a  power  study  would  need  to  be  performed 
to  determine  when  sample  sizes  are  large  enough. 
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